home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / interpolation / poly.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-12-03  |  4.2 KB  |  174 lines
  1. /* interpolation/interp_poly.c
  2.  * 
  3.  * Copyright (C) 2001 DAN, HO-JIN
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_errno.h>
  23. #include <gsl/gsl_poly.h>
  24. #include <gsl/gsl_interp.h>
  25.  
  26. typedef struct
  27. {
  28.   double *d;
  29.   double *coeff;
  30.   double *work;
  31. }
  32. polynomial_state_t;
  33.  
  34. static void *
  35. polynomial_alloc (size_t size)
  36. {
  37.   polynomial_state_t *state =
  38.     (polynomial_state_t *) malloc (sizeof (polynomial_state_t));
  39.  
  40.   if (state == 0)
  41.     {
  42.       GSL_ERROR_NULL ("failed to allocate space for polynomial state",
  43.               GSL_ENOMEM);
  44.     }
  45.  
  46.   state->d = (double *) malloc (sizeof (double) * size);
  47.  
  48.   if (state->d == 0)
  49.     {
  50.       free (state);
  51.       GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
  52.     }
  53.  
  54.   state->coeff = (double *) malloc (sizeof (double) * size);
  55.  
  56.   if (state->coeff == 0)
  57.     {
  58.       free (state->d);
  59.       free (state);
  60.       GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
  61.     }
  62.  
  63.   state->work = (double *) malloc (sizeof (double) * size);
  64.  
  65.   if (state->work == 0)
  66.     {
  67.       free (state->coeff);
  68.       free (state->d);
  69.       free (state);
  70.       GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
  71.     }
  72.  
  73.   return state;
  74. }
  75.  
  76. static int
  77. polynomial_init (void *vstate,
  78.          const double xa[], const double ya[], size_t size)
  79. {
  80.   polynomial_state_t *state = (polynomial_state_t *) vstate;
  81.  
  82.   int status = gsl_poly_dd_init (state->d, xa, ya, size);
  83.  
  84.   return status;
  85. }
  86.  
  87. static int
  88. polynomial_eval (const void *vstate,
  89.          const double xa[], const double ya[], size_t size, double x,
  90.          gsl_interp_accel * acc, double *y)
  91. {
  92.   polynomial_state_t *state = (polynomial_state_t *) vstate;
  93.  
  94.   *y = gsl_poly_dd_eval (state->d, xa, size, x);
  95.  
  96.   return GSL_SUCCESS;
  97. }
  98.  
  99.  
  100. static int
  101. polynomial_deriv (const void *vstate,
  102.           const double xa[], const double ya[], size_t size, double x,
  103.           gsl_interp_accel * acc, double *y)
  104. {
  105.   polynomial_state_t *state = (polynomial_state_t *) vstate;
  106.  
  107.   gsl_poly_dd_taylor (state->coeff, x, state->d, xa, size, state->work);
  108.  
  109.   *y = state->coeff[1];
  110.  
  111.   return GSL_SUCCESS;
  112. }
  113.  
  114. static int
  115. polynomial_deriv2 (const void *vstate,
  116.            const double xa[], const double ya[], size_t size,
  117.            double x, gsl_interp_accel * acc, double *y)
  118. {
  119.   polynomial_state_t *state = (polynomial_state_t *) vstate;
  120.  
  121.   gsl_poly_dd_taylor (state->coeff, x, state->d, xa, size, state->work);
  122.  
  123.   *y = 2.0 * state->coeff[2];
  124.  
  125.   return GSL_SUCCESS;
  126. }
  127.  
  128. static int
  129. polynomial_integ (const void *vstate, const double xa[], const double ya[],
  130.           size_t size, gsl_interp_accel * acc, double a, double b,
  131.           double *result)
  132. {
  133.   polynomial_state_t *state = (polynomial_state_t *) vstate;
  134.   size_t i;
  135.   double sum;
  136.  
  137.   gsl_poly_dd_taylor (state->coeff, 0.0, state->d, xa, size, state->work);
  138.  
  139.   sum = state->coeff[0] * (b - a);
  140.  
  141.   for (i = 1; i < size; i++)
  142.     {
  143.       sum += state->coeff[i] * (pow (b, i + 1) - pow (a, i + 1)) / (i + 1.0);
  144.     }
  145.  
  146.   *result = sum;
  147.  
  148.   return GSL_SUCCESS;
  149. }
  150.  
  151. static void
  152. polynomial_free (void *vstate)
  153. {
  154.   polynomial_state_t *state = (polynomial_state_t *) vstate;
  155.   free (state->d);
  156.   free (state->coeff);
  157.   free (state->work);
  158.   free (state);
  159. }
  160.  
  161. static const gsl_interp_type polynomial_type = {
  162.   "polynomial",
  163.   3,
  164.   &polynomial_alloc,
  165.   &polynomial_init,
  166.   &polynomial_eval,
  167.   &polynomial_deriv,
  168.   &polynomial_deriv2,
  169.   &polynomial_integ,
  170.   &polynomial_free,
  171. };
  172.  
  173. const gsl_interp_type *gsl_interp_polynomial = &polynomial_type;
  174.